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We investigate the energy transport in one-dimensional disordered granular solids by extensive 
numerical simulations. In particular, we consider the case of a polydisperse granular chain composed 
of spherical beads of the same material and with radii taken from a random distribution. We start 
by examining the linear case, in which it is known that the energy transport strongly depends on the 
type of initial conditions. Thus, we consider two sets of initial conditions: i) an initial displacement 
and ii) an initial momentum excitation of a single bead. After establishing the regime of sufficiently 
strong disorder, we focus our studies on the role of nonlinearity for both sets of initial conditions. 

By increasing the initial excitation amplitudes we are able to identify three distinct dynamical 
regimes with different energy transport properties: a near linear, a weakly nonlinear and a highly 
nonlinear regime. Although energy spreading is found to be increasing for higher nonlinearities, in 
the weakly nonlinear regime no clear asymptotic behavior of the spreading is found. In this regime, 
we additionally find that energy, initially trapped in a localized region, can be eventually detrapped 
and this has a direct influence on the fluctuations of the energy spreading. We also demonstrate that 
in the highly nonlinear regime, the differences in energy transport between the two sets of initial 
conditions vanish. Actually, in this regime the energy is almost ballistically transported through 
shock-like excitations. 


I. INTRODUCTION 

Wave scattering and energy transport in disordered 
media have been for long time a matter of great research 
interest [1]. The experimental observation of Anderson 
localization in different systems such as optical [2] , ultra¬ 
cold atomic gases [3] and elastic networks [4] , has renewed 
the research in this direction. In addition, recent stud¬ 
ies on wave scattering in random media have lead to a 
plethora of applications in imaging, random lasing and 
solar energy (see for example [5] and references therein). 

The key phenomenon employing in these studies is the 
spatial wave localization due to disorder, which is a linear 
effect relying on keeping phase coherence of participat¬ 
ing waves [6]. However, wave localization can also emerge 
due to nonlinearity, as it was first shown in the studies 
of Fermi-Pasta-Ulam (FPU) [7], and may lead to energy 
localization and propagation through the formation of 
localized solutions (solitons, breathers etc.) in different 
lattice models [8] . The interplay of these two localization 
mechanisms, nonlinearity and disorder, has been studied 
extensively in the recent years [9H26]. In most of these 
studies, an initially localized wavepacket was shown to 
lead to delocalization and a sub-diffusive spreading of 
the energy, for sufficiently large nonlinearities. The most 
common models that have been studied are the Klein- 
Gordon (KG) model as well as the discrete nonlinear 
Schrodinger (DNLS), where especially the latter has at¬ 
tracted much attention due to its application to various 
optical structures and devices. Experimental studies on 
optical structures show that nonlinearity can either en¬ 
hance localization (for focusing nonlinearity) or induce 
delocalization (for defocusing nonlinearity) [27l |28] . 

Granular solids, namely densely packed arrays of 
macroscopic particles which appear naturally disordered, 


are a promising testbed for studying the interplay of dis¬ 
order and nonlinearity. The latter originates from the in¬ 
terparticle Hertzian contacts [29]. An especially appeal¬ 
ing characteristic of these media is their tunable dynami¬ 
cal response ranging from near linear to highly nonlinear, 
by changing the ratio of static to dynamic interparticle 
displacements. Fabricated granular solids have allowed 
the exploration of a plethora of fundamental phenomena, 
including solitary waves with a highly localized waveform 
in the case of uncompressed crystal, discrete breathers 
and others [SOHSZ]- They have been also applied in var¬ 
ious engineering devices, including among others shock 
and energy absorbing layers [38ti4Q] . acoustic lenses m, 
and acoustic diodes [42] . 

For sufficiently weak excitations and in the presence of 
precompression, the one-dimensional disordered granular 
solid, also called granular chain, can be approximated by 
a disorder harmonic lattice which has some interesting 
transport properties. In particular, it has been shown 
that different initial conditions - initial displacement or 
momentum excitations - of a single particle can lead to: 
sub-diffusive (displacement) or super-diffusive (momen¬ 
tum) energy transport [43], as well as to analytical de¬ 
scribed asymptotic energy profiles m- On the other 
hand, for sufficiently strong excitations or in the absence 
of precompression, the granular chain exhibits two differ¬ 
ent types of nonlinearity: (i) a power nonlinearity stem¬ 
ming from the Hertzian contacts and (ii) a non-smooth 
nonlinearity, which is triggered whenever two beads of 
the chain lose contact (gap opening). The latter is 
present to a broad class of fragile mechanical systems 
that loose rigidity upon lowering the external pressure 
towards zero, such as weakly connected polymers [53| 
and network glasses [45]. It is also present in cracked 
solids [46] . 
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Recently, studies of one-dimensional disordered gran¬ 
ular chains have been reported [UHH]. In the absence 
of precompression, where the non-smooth nonlinearity is 
present, it was shown that if a solitary wave is formed, it 
features an exponential decay which strongly depends on 
the degree of randomness [48l [49] . Similar results were 
also reported in a two-dimensional granular solid [50| 
where the decay of the amplitude of the wave front was 
described using an analogy between disorder and vis¬ 
coelastic dissipation. On the other hand, in the pres¬ 
ence of precompression, the power nonlinearity stemming 
from the Hertzian contacts leads to a FPU like dynamics, 
which have been studied theoretically in the presence of 
disorder [24H26] . However, in the case of granular chains, 
this dynamics can be strongly modified by the presence 
of the opening of gaps. Thus, the interplay of these two 
nonlinear mechanisms is of particular interest and it can 
drastically change the transport properties. 


Only recently, a study about one-dimensional and pre¬ 
compressed random dimer granular chains [51] has re¬ 
ported some features of the energy transport. In this 
work, the authors compare wave dynamics in chains 
with three different types of disorder: an uncorrelated 
(Anderson-like) disorder and two types of correlated dis¬ 
order. For the Anderson-like uncorrelated disorder, they 
found a transition from subdiffusive to superdiffusive dy¬ 
namics depending on the amount of precompression in 
the chain. In the present work, we consider a differ¬ 
ent kind of uncorrelated disorder (i.e. polydispersity 
through disorder in the bead radius) and we study both 
displacement and momentum initial excitations, empha¬ 
sizing their differences and similarities [52]. In particular, 
we consider polydisperse disordered granular chains com¬ 
posed of spherical beads of the same material and with 
radii taken from a random distribution. Our motivation 
is the fact that most of the granular materials occurring 
in nature and industrial application are composed of a 
broad range of particle sizes [53] . By considering a single 
central bead excitation, we study the transport of energy 
in these disordered granular chains. 


In Sec. [H] we present the equations of motion in a nor¬ 
malized form, we define the conserved energy of the sys¬ 
tem and also describe the parameters used for the char¬ 
acterization of the energy transport. Results for the lin¬ 
ear case are shown in Sec. im where the influence of the 
strength of the disorder on the dynamics is studied. In 
Sec. lEl we show the main results of this work for the 
case of an initial displacement excitation, and explic¬ 
itly identify three different regions of energy transport, 
a near linear one, a weakly nonlinear and a highly non¬ 
linear regime. The energy transport for the case of an 
initial momentum excitation is discussed in Sec. |V| pre¬ 
senting differences and similarities with the initial dis¬ 
placement excitation. Results concerning the asymptotic 
profile of the energy for the two types of initial conditions 
are shown in Sec. VI and finally in Sec. VH we conclude 
our results. 
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FIG. 1. A sketch of a granular chain with beads of random ra¬ 
dius. Un denotes the displacement of each bead from its equi¬ 
librium position, while Sn is the overlap between two spherical 
beads due to the precompression force F. 


II. DISORDERED GRANULAR CHAIN 

We consider a one-dimensional chain consisting of V + 
2 spherical beads, with masses (n = 0,1, 2 ... V + 1) 
and Hertzian contacts as shown in Fig. We consider 
fixed boundary conditions for the first and last spherical 
bead, namely uq = un-\-i = 0, where Un is the displace¬ 
ment of each bead from its equilibrium position. Then, 
the system is described by the following set of differential 
equations: 


3 3 

miili = Ai[6i — ui]^ — A2[S2 ui — + , 

3 3 

'^nhn — [^n 1 '^n]-|- ^n+1 [^n+1 '^n+l]-|-5 

3 3 

tunUn = An[Sn + un-1 — un]^ — An^i[Sn^i + un]^, 

(1) 

where A^ is the contact coefficient between beads n — 1, n 
and Sn is the relative static overlap due to a precompres¬ 
sion force F acting on the two boundaries. The dots (') 
denote differentiation with respect to time. The coeffi¬ 
cient An for spherical beads of the same material, is given 

by An = (2/3)£^(Rn-lRn)/(En-l + En)/(1 “ 122], 

where F and n are the elastic modulus and Poisson’s ra¬ 
tio respectively, while Rn is the radius of the nth bead. 
The static overlap Sn is given by Sn = (F/An)^^^ [29] . 
The [ ]+ sign in Eq. 0 denotes the following: when 

the expression inside the square brackets becomes nega¬ 
tive (i.e the beads are not in contact) this term becomes 
zero. In fact this happens when the relative displace¬ 
ment between two beads becomes larger than their over¬ 
lap Un-i —Un>Sn, that is there is a gap between them, 
and their relative force vanishes. 

Below we will work in dimensionless units, however 
for clarity we note that we use a reference radius of 
R = 0.01 m, a static force of F = 1 N and stainless 
steel spherical beads (316 type), the elastic modulus of 
which is £ = 193 GPa while the Poisson ratio is u = 0.3. 
Relevant experiments with granular chains contained few 
defects can be found in [54]. In the following we will 
consider a disordered setup where the radii Rn of the 
different beads will be taken as a random variable, with 
values taken from a uniform distribution within the range 
Rn G [R, aR] , where the parameter a > 1 describes the 
disorder strength. Consequently, the mean value R of 
the bead radius is R = {a -\- l)R/2. In order to make 
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our equations dimensionless we implement the following 
transformations for time, distance, mass and stiffness re¬ 
spectively: 

t UJt^ 6n {Un . . 

^ \^) 

TTin ^ ^ 


where all the quantities with tilde, are calculated at R. 
The frequency uJc = is the cutoff fre¬ 

quency corresponding to the linear case of a chain with 
spherical beads of radius R. The normalization is such 
that in the case of no disorder (a = 1) the normalized 
cutoff frequency is cJc = 1- The energy of the system is 
given by the following expression: 

N ^ / 2 \ 

n=l n=l \ 'Ti J 

where and Pn = rUnUn are the energy and momentum 
of the nth bead respectively. The potential Vn for each 
spherical bead is defined as Vn = \y{un) + y(n^+i)]/2 
where: 

V{Un) = ‘^An[5n + Mn-l “ Mn]+^ “ 

^n^n iVn—l "^n)* 

To study the energy transport in this one-dimensional 
system we focus on the time evolution of the second mo¬ 
ment of the energy distribution [26] defined as: 

n,, = (5) 

E 

where no = N/2 corresponds to the central bead of the 
chain and N = 2 x 10^ to the total number of the spherical 
beads. In the last expression, hn = EnjE denotes the 
portion of the total energy E acquired by the nth bead. 
Another useful quantity that characterizes the system is 
the participation number: 

P=l/Y,hl (6) 

n 

which measures the number of excited beads that signif¬ 
icantly contribute in the energy distribution. It takes the 
value P = 1 if all the energy is concentrated in one bead 
while it becomes P = N in the case of energy equipar- 
tition. In our study, we investigate the energy transport 
under two different sets of initial conditions: 


• WiV/2(0) = u, U„^JV/2(0) = 0, Un{0) = 0, 

• Mn(0) = 0, UiV/2(0) = U, Un^N/2{0) = 0, 

corresponding respectively to an initial displacement and 
an initial momentum excitation of the central bead. We 
present results obtained by averaging over 200 disorder 
realizations. Throughout the text the average value over 
disorder realizations of a quantity x, is denoted as {x). 
Simulations are carried out using the SABA2C symplec- 
tic integrator which allows as to keep the relative energy 
error at the order of 10“^ [ming. We also note that in 
our simulations energy never reaches the boundaries of 
the chain. 



ol_ ^^^ ^ 
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FIG. 2. (Color online) Harmonic chain: (a) Time evolu¬ 
tion of (log7712(t)) for different values of a. The shaded 
area corresponds to the standard mean deviation of the mea¬ 
sured mean value. Curves from top to bottom correspond 
to o = 1.1,1.5, 2,3 respectively, (b) The time derivative /3 
of (logiom 2 (t)) given by Eq. as a function of log^Q t for 
a — 1.1,1.5, 2,3 from left to right. The horizontal dashed 
lines indicate the value /3 = 0.5. (c) Time evolution of the 
mean participation number (P) for the different values of the 
disorder parameter a. Curves from top to bottom correspond 
to a — 1.1,1.5, 2, 3 respectively, (d) The asymptotic value of 
the mean participation number (P)* as a function of a. 


III. HARMONIC CHAIN 

For sufficiently small displacements, i.e. Un-i —Un^ 
Sn, the system of Eqs. 0 can be approximated by the 
following linear system: 

'^n'dn — Eni'^n—l A^n+l('^n "^ 71 + 1)5 (^) 

where Kn = {3/2)AnSn^‘^ is the linear coupling constant. 
In the absence of disorder, a = 1, the energy spread¬ 
ing is ballistic and the second moment grows in time as 
m 2 {t) ^ while (P) diverges for both displacement and 
momentum initial excitations. On the other hand, for the 
case of randomly chosen radii, Eq. 0 has the form of a 
disordered harmonic chain. This system has already been 
studied in several works [26l |43l [56] with either a mass 
disorder or a disorder in the coupling constants An- For 
the granular chain considered here, having beads of the 
same material but of different radius, both the masses 
mn and the coupling constants Kn are random variables 
(both depend on the radius of the beads). Since masses 
depend on the radius as m oc P^, while the linear cou¬ 
plings as A oc P^/^,we expect that the disorder effect is 
stronger due to the masses. 

In order to investigate the importance of the disorder 
parameter a on the system’s behavior, we numerically 
integrate Eqs. 0 for different values of a and the results 
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are presented in Fig. In particular, in Fig. (a), we 
show the time evolution of the average logarithm of the 
second moment (log^^Q m 2 (t)) with respect to the loga¬ 
rithm of time. Furthermore in the panels of Fig. i(b), 
we show the time evolution of yd, the time derivative of 
(logic m 2 (i)) given by 


^(logiomz) 
d logic i 


(8) 


The derivative is calculated numerically as follows: we 
first smooth the values of (logio(^ 2 )) by using a lo¬ 
cally weighted regression algorithm [57], and then we 
apply an 8^^ order central finite difference scheme to 
compute the derivative. In [26l [43], it was shown that 
for an initial displacement the energy transport is sub- 
diffusive. In particular, the second moment grows in time 
as m 2 (t) ^ with an asymptotic value for the exponent 
yd(t ^ oo) ^ 0.5. From the left panel of Fig. (b), 
its readily seen that for a = 1.1 the parameter yd(t) ini¬ 
tially acquires a value close to yd = 2 indicating a ballis¬ 
tic spreading of energy, but eventually it drops to smaller 
values implying a slower spreading. However, one can not 
induce a clear sub-diffusive behavior for the time scales 
of our simulations. The second and third panels of Fig.|^ 
(b) correspond to stronger disorder {a = 1.5 and a = 2 
respectively) where the tendency of /3 to asymptotically 
reach the value of yd = 0.5 is evident, although some 
larger fluctuations are present. In the rightmost panel of 
Fig. i (b), we plot the case of a = 3, we see that the 
value of yd saturates to yd = 0.5 very fast. However, we 
also notice that for this value of a, the fluctuations are 
even larger and this is also depicted by the very large 
standard deviation in the mean value of {log^Q m 2 {t)), 
indicated by the shaded area around the lowest curve of 
Fig. I (a). Therefore, a much larger number of realiza¬ 
tions is needed for a better analysis of the a = 3 case. 

The respective mean participation number (P) for the 
different values of the disorder parameter a is shown in 
Fig. i (c). For a = 1.1, the mean participation num¬ 
ber does not saturate, at least on the times of our sim¬ 
ulations. On the other hand, for stronger disorder e.g. 
(a = 2, it acquires an asymptotic value depending on the 
disorder parameter. The dependence of the numerical es¬ 
timation of the asymptotic value of (P(t ^ oo)) = (P)* 
for different values of a is shown in Fig. (d). These 
estimations are obtained using the results of Fig. (c) 
as the mean value of (P) at the second half of the last 
decade of the simulation i.e. for 5 x 10^ < t < 10^. Thus 
we may identify: a weak disorder regime {a < 1.5) where 
(P)* has a value larger than 10 and approaches the total 
number of particles N in the case of no disorder {a = 1), 
and a strong disorder regime {a > 1.5) where it satu¬ 
rates to values of about 10 beads or less. For this reason 
and due to the fact that for values a > 3 the smoothing 
of large fluctuations of the computed quantities would 
imply many more disorder realizations, we restrict our 
analysis to the disorder regime with a = 2. 

We note that similar results to the ones of Fig. were 
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FIG. 3. Time evolution of the logarithm of the energy En 
of a portion of the chain near the initially excited bead at 
n = N/2 = 10^. Different panels (from top to bottom) corre¬ 
spond to different values of the initial displacement u, and in 
particular to = 0.01, 1, and 10 respectively. The coloring of 
each lattice site, according to the color scale shown on top of 
the panels, denotes the log^Q (Pn) value of the corresponding 
bead. The insets on the right of each panel show the corre¬ 
sponding energy profile at a late time instant, t = 2 x 10^. 


obtained for the case of initial momentum excitation. In 
particular, we recovered that the asymptotic value of [3 
is 1.5, in the case of a = 2. In the rest of this work we 
systematically investigate the effect of the nonlinearity 
to our system. 


IV. DISPLACEMENT EXCITATION 

In this Section we present our findings for the trans¬ 
port of energy induced by an initial displacement of the 
central spherical bead in the chain. Typical results of the 
dynamics observed for different amplitudes of the initial 
displacements, are shown in Fig. The top and middle 
panels illustrate clearly that a localized state is formed 
near the initially excited bead, while there are two prop¬ 
agating fronts traveling towards the edges of the chain. 
These results show that for an increasing amplitude of 
the initial excitation, these fronts are found to be prop¬ 
agating faster. In the bottom panel, which corresponds 
to 14 = 10, the behavior is significantly altered, since now 
the energy looks to be more equally distributed through 
the lattice. This fact is better illustrated by looking at 
the insets (in the right of each panel), where we show 
the corresponding energy profile at a late time instant 
t = 2 X 10^ for each case. For the cases of the top and 
middle panels, it is readily seen that the energy around 
the initially excited bead is almost five orders of magni¬ 
tude larger compared to the energy of more distant beads. 
On the other hand, in the bottom panel the differences 
of energy between the central region and the rest of the 
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FIG. 4. (Color online) Time evolution of: (a) the averaged 
logarithm of the second moment (log^o ^ 2 ) for different initial 
displacements (b) the time derivative /3(t) of the respective 
curves of panel (a), (c) the mean value of the participation 
number (P). 


eter [3 for different values of the initial displacement excita¬ 
tion u = 0.2, 0.4, 0, 6 (from left to right), in the regime where 
energy transport crosses from sub- to super-diffusive behav¬ 
ior. The horizontal dashed lines indicate the value /3 = 1. 
(b) Time evolution of the mean participation number for dis¬ 
placement excitations corresponding to the top panels. 


chain are significantly smaller. 

In what follows we discuss in more detail the outcomes 
of extensive numerical simulations for several values of 
the nonlinearity strength. The main results are summa¬ 
rized in Fig. In particular, in Fig. 0 (a) we plot the 
time evolution of (log^^o ^ 2 ), in Fig.|^(l^we plot the time 
evolution of its time derivative and finally in Fig.[^(c) we 
show the mean participation number (P) as a function 
of time. 


A. Near linear regime 

Let us first note that for values of the initial excita¬ 
tion u < 0.1 [see e.g. the blue (lowest) line in Fig. El(b) 
which corresponds to u = 0.01], we observe that p ^ 0.5 
at about t ~ 10^. Additionally the mean participation 
number for this case [Fig. (c)] is found to practically 
saturate to the same value as in the harmonic chain. 
Thus, we conclude that for values oi u < 0.1 the non¬ 
linear model is characterized by a near linear regime, 
having qualitatively and quantitatively the same energy 
transport properties as its linear counterpart, at least for 
the time scales of our simulations. 


B. Weakly nonlinear regime 

From the results of Fig. it is readily seen that already 
for an initial condition of rt = 0.1, the value of (3 deviates 
from its behavior in the linear case, becoming larger than 
0.5. However a very abrupt change in the behavior of p is 
observed at the value u = 1 and the system undergoes a 
transition from sub- to super-diffusion. Notice also that 
although for u = 0.1 the mean participation number in 


Fig. 0 (c), saturates to a value of (P) ~ 10, for the value 
14 = 1 it is found to continuously increase. These two 
results indicate the existence of a regime of intermediate 
dynamics for values between i4 = 0.1toi4 = l. 

In order to further understand the dynamics in this 
regime, we compute the parameter ^5, as well as the mean 
participation number (P), for some intermediate values 
of initial displacements i.e. u = 0.2,0.4,0.6. The ob¬ 
tained results, which are plotted in Fig. clearly show 
that a transition from sub- to super-diffusion is carried 
out in this regime. In Fig. (a), the parameter /3 ex¬ 
hibits many fluctuations and shows no evident tendency 
to saturate into a constant value until the end of our nu¬ 
merical simulations. In all cases shown in Fig. Ha), /3 
initially approaches the diffusive value P = I but later 
on starts to decrease. Furthermore at a time interval be¬ 
tween t ~ 10^ and t ~ 10^ it saturates to an almost 
constant value somewhat below /3 = 1, but eventually 
the dynamics changes and P starts to increase again be¬ 
coming larger than 1. This behavior creates a character¬ 
istic local minimum of P{t) for all studied cases shown 
in Fig. This result is in accordance with the recent 
study of [26], where it was found that in the FPU prob¬ 
lem, until the end of the studied times, there was no clear 
asymptotic value for the exponent of m 2 . 

It is also relevant to discuss the behavior of the 
mean participation number (P) in this weakly nonlinear 
regime. As it is shown in Fig. (b) after an initial in¬ 
crease of (P) its value remains practically constant with 
a value of around (P) ^ 10, between t ^ 10^ and t ^ 10^, 
but after this time interval (P) increases exhibiting a di¬ 
verging trend. This transition is attributed to nonlinear¬ 
ity, since it is never observed in the near linear regime or 
in the exact linear case. To further understand the origin 
of this behavior, we plot in Fig. |^the normalized mean 
energy of the central bead {hjsf/ 2 ) as a function of time. 
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FIG. 6. (Color online) The normalized mean energy {hN/ 2 ) 
of the central bead at n = A/'/2 as a function of time, for 
different initial displacements; from bottom to top, u = 
0 . 1 , 0 . 6 , 1 , 10 , 100 . 


For rt = 0.1 after a small transient time of t ~ 10, the 
central bead retains a large amount of the total energy 
of the system, keeping it up to the end of our simula¬ 
tions at t = 10^. This behavior is understood by the 
fact that most of the energy is concentrated around the 
initially excited bead, due to the presence of the disorder- 
induced localized modes. These modes remain localized 
throughout the simulation. However, for larger values 
of the initial excitation, i.e. 0.1 < rt < 1, we observe that 
although for a large time interval (10 < t < 10^) most 
of the energy is trapped around the central bead, after 
sufficient time the energy of this bead starts to decrease. 
This signals the detrapping of energy from this bead, and 
its release to the rest of the chain. 


C. Highly nonlinear regime 

For even larger values of the initial condition, u > 
i.e. for larger nonlinearities, we observe in Fig. that 
the exponent /3 saturates to an almost constant value at 
about t ~ 10^ describing a super-diffusive regime. These 
values are larger with respect to the values of (3 seen in 
the weakly nonlinear regime. It is worth noting that in 
this regime the fluctuations in the values of (3 are much 
smaller than in the previous regime. From Fig.|^ we also 
conclude that the normalized mean energy of the central 
bead (/iAr/ 2 ) foi* = 10? 100 continuously decreases as 
a function of time, in contrast with the weakly nonlin¬ 
ear regime. To further investigate the dynamics in this 
regime we evaluate the probability of gap openings be¬ 
tween beads as obtained by counting the number of gaps 
in each site for all the 200 disorder realizations, and plot 
in Fig.j^the obtained average value. In great contrast to 
the weakly nonlinear regime where no gaps appear, here 
we find that not only there are always many gaps around 
the initially excited bead, but also these gaps propagate 
in the system. This observation strongly suggests that 
for such large initial excitations, a new dynamical regime 
is present, where the dynamics is not governed by the 
FPU like nonlinearity but by the nonsmooth nonlinear¬ 
ity of the opening of gaps. We call this regime highly 



FIG. 7. The probability of a gap opening as a function of time 
for two different initial displacements in the highly nonlinear 
regime: u = 10 (upper panel) and u = 100 (lower panel). 

We focus on the dynamics around the central bead N/2. The 
black color in the colormap corresponds to probability 1 while 
the white to 0. 



FIG. 8. (Golor online) The velocity prohles iin for hve time 
instants, normalized to its maximum value for different initial 
conditions. Left: the case of T = 0 with u = 1, middle and 
right: the case of F ^ 0 with = 100 and u = 10 respectively. 


nonlinear. 

A particular limit of this highly nonlinear regime cor¬ 
responds to the case of F = 0 which results to = 0 in 
Eq. 0 - This is also called sonic vacuum 1201 due to the 
fact that the system does not support the propagation of 
linear waves. This regime has been studied extensively 
for the case of no disorder, namely when a = 1. It is 
known that a solitary solution exists in this limit with a 
highly localized waveform m- For the case of a binary 
system with a disordered distribution between beads of 
two different masses, similar solitary waves of decreasing 
amplitude were found in the weak disorder limit, while 
in the strong disorder case a delocalized wave was ob¬ 
served [ 49 ]. Similar results were obtained for the case of 
two-dimensional granular solids, with an initial excita¬ 
tion only in one direction: weak disorder induces an expo¬ 
nentially decreasing solitary wave which eventually gives 
its place to a delocalized shock-like profile, while strong 
disorder only exhibits the shock-like structure [50l [58] . 
The latter works also showed that at the position Uf of 

the front of the shock-like structure, the velocity scales 
-1/2 

as Uf (X 

Let us now explore the transition from the highly non¬ 
linear regime to the singular case of F = 0. First note 
that, as it is shown in Fig. (b), for F = 0 the exponent 
(3 reaches the value ~ 1.8, which is the maximum value. 
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FIG. 9. (Color online) Same as in Fig. but for different 
initial momentum excitations un/ 2 { 0 ) = ii. 



FIG. 10. [Color online] Same as in Fig. but for different 
initial momentum excitations '?iiv/2(0) = ii. The top panels 
correspond to = 0.1, 0.2, 0.4 from left to right. 


while the mean participation number is qualitatively sim¬ 
ilar for all the cases in the highly nonlinear regime. Fur¬ 
thermore in Fig. we plot the velocity profiles at five 
time instants for the case of F = 0 with u = and for 
the case of F 7 ^ 0 with u = 10 and u = 100. For F = 0, 
in agreement with [50l [58] , we observe the formation of 
a shock-like structure, see left case in Fig. p] In fact, 
by plotting a fitting curve of the form , shown as 

an envelope on-top of the velocity profiles [solid (green) 
line], it is clearly seen that the propagating front follows 
this — 1/2 power law trend. 

More importantly we find that even in the case of a 
finite precompression force, F 7 ^ 0, a similar structure 
can be formed and is propagating with the same power 
law decay of its amplitude, see middle case in Fig. 
However, for the case of u = 10 although the peak of 
the propagating form seems to follow a similar trend, 
the observed structure is different; instead of a shock¬ 
like profile it appears to be bell-shaped and it consists of 
a smaller number of beads (~ 100 beads), see right case 
in Fig. 


V. MOMENTUM EXCITATION 

As we mentioned in the introduction, the energy trans¬ 
port in disordered linear chains, strongly depends on the 
initial condition [26l [43] . Thus in order to complete the 
dynamical study of our model we also investigate the case 
of an initial momentum excitation of the central spher¬ 
ical bead by increasing the initial velocity Un/ 2 {^) = u. 
The results are summarized in Fig. in and Fig. 

In the linear case the energy spreading for an initial 
momentum excitation is known to be super-diffusive with 
(3 = 1.5 [43]. By numerical simulations of the normalized 
nonlinear equations of motion (^, we found that the en¬ 
ergy spreading remains super-diffusive with (3 = 1.5, for 
all initial velocities -h < 10 [see Fig. |^(b), Fig.[T^(a)]. 


This is in contrast with respect to the case of the initial 
displacement excitation, where the relevant values of (3 
exhibit large variations. For -h > 10 the exponent (3 grows 
and eventually reaches the maximum value of f3 ^ 1 . 8 , 
which is also the corresponding value of (3 of the limit¬ 
ing case of F = 0. We again identify this regime as the 
highly nonlinear regime. 

On the other hand, the mean participation number (F) 
exhibits similar behavior with the case of initial displace¬ 
ment excitations. For -u < 0.1 it saturates to a constant 
value of about 20 beads, while for > 0.1 it continuously 
increases. Looking closer to (F), shown in Fig. [lQ|[b) 
we note that the mean participation number reaches the 
asymptotic constant value of about 20 beads for -h = 0 . 1 , 
i.e., the near linear regime. However, for it = 0.2, 0.4 al¬ 
though it seems to saturate to a constant value for long 
time, finally it starts to deviate and to increase continu¬ 
ously. This again suggests that there is an intermediate 
regime for 0.1 < ft < 1 in which energy detrapping is 
observed. This is the weakly nonlinear regime. 

To conclude, compared to the case of displacement ini¬ 
tial excitations, there is a difference in the energy trans¬ 
port properties in the intermediate regime that we call 
weakly nonlinear. Although the mean participation num¬ 
ber (F) shows the same behavior, in contrast to the dis¬ 
placement excitations, for momentum initial excitations 
both the near linear and weakly nonlinear regimes are 
characterized by an asymptotic value of the parameter 
(3 = 1.5, namely the same as to the linear case. 

Another interesting point to be mentioned is that al¬ 
though in the linear case the energy transfer is sub- 
diffusive or super-diffusive for an initial displacement 
or an initial momentum excitation respectively, in the 
highly nonlinear regime and more profoundly in the lim¬ 
iting case of F = 0 both excitations result in the same 
behavior of energy transport. In fact in this limit, not 
only the asymptotic value of /3 (as it was also found in 
Ref. [52] for a disordered dimer), but also the dynamics 
of the derivative P of (m 2 ) are very similar. 






















VI. ASYMPTOTIC PROFILES 


According to the recent work of Ref. m, the asymp¬ 
totic dependence of the energy moments (not only m 2 ) 
can be characterized by the asymptotic energy profile of 
the lattice, at least in the linear case. In particular, it was 
shown in that work, that the energy profile far away from 
the central excitation has the form of (hn) ^ \n — N/2\~'^. 
The exponent 77 was found to be 77 = 5/2 and 77 = 3/2 for 
displacement and momentum initial excitations respec¬ 
tively. In Fig. 11 we show for both displacement (top 


panels) and momentum initial excitations (bottom pan¬ 
els) three different instants of the energy profile at suf¬ 
ficiently large times, covering all the dynamical regimes 
from near linear (left column) to highly nonlinear (right 
column). In the left panels (near linear regime) it is read¬ 
ily seen that, in accordance with the predictions for the 
linear problem, the three profiles overlap indicating the 
fact that the energy distribution has reached an asymp¬ 
totic profile which does not change for later times. For 
comparison, we also plot curves (dashed lines) with a 
slope of 5/2 (top) and 3/2 (bottom). Note also that 
the energy distribution near the central bead (i.e. for 
n ~ N/2 in the figure) is the same for the three dif¬ 
ferent profiles. This indicates that the energy of these 
sites also does not change in time, and the mode around 
the center remains localized. 

The top (bottom) second column panel of Fig. de¬ 
picts the energy profile for an initial displacement (mo¬ 
mentum) excitation with u = 1 {ii = 1). In this case it is 
readily seen that for the initial displacement excitation 
(top) not only the three profiles do not overlap, but also 
the slopes are not exactly 5/2. This is another indication 
for the appearance of the weakly nonlinear regime which 
exhibits nontrivial dynamics. Note that for the momen¬ 
tum excitation (bottom) the profiles do overlap. This is 
in accordance with the discussion in section V about the 
asymptotic value of (3 which is the same for both the near 
linear and the weakly nonlinear regime. Additionally, for 
both cases of initial conditions, we observe a small but 
non negligible deviation of the energy around the center, 
which confirms the fact that the localized mode around 
the center starts to loose its energy. 

For larger values of the initial condition as is shown 
in the top (bottom) third column panel of Fig. 11, the 
profile of the energy is substantially different. It is char¬ 
acterized by two different regimes: a weakly localized 
part with almost 10 ^ beads with similar amount of en¬ 
ergy, forming an almost straight horizontal line and a 
decaying tail. However, the weakly localized portion of 
the energy is spreading and loses its energy, as for larger 
times the almost straight horizontal part of the profile 
becomes longer having smaller energy values. We also 
note that once again the slope of the energy profiles for 
initial displacements (top panel) is not 5/2 as time in¬ 
creases but interestingly enough it reaches a value which 
is closer to 3/2 (see gray solid line). We remind that 
for this initial displacement excitation with tx = 10, /3 


asymptotically reaches the value of 1.5 (see Fig.[^ which 
is also the asymptotic value of [3 of an initial momentum 
excitation of the linear problem. 

Finally for the particular case of F = 0, as shown in the 
top (bottom) fourth column panel of Fig. 11, the asymp¬ 


totic energy profile is similar to a ballistic propagation 
where the energy differences between excited beads de¬ 
crease drastically but in this case the propagating front 
does not exhibit a sharp profile. In fact this front has a 
very large “tail” which is due to the shock-like structure 
that was mentioned in section C. 


VII. CONCLUSIONS 

In this work we numerically investigated the energy 
transport in a one-dimensional granular solid composed 
of spherical beads of randomly distributed radii which 
interact via Hertzian forces. We studied the dynamics 
by using two different localized initial conditions i.e. ini¬ 
tial displacement and initial velocity excitations of the 
central bead of the chain and by increasing the ampli¬ 
tude of these excitations. We were able to identify three 
different dynamical spreading regimes with distinct char¬ 
acteristics: the near linear, the weakly nonlinear and the 
highly nonlinear. 

In the near linear regime, part of the initial energy 
remains localized around the central excited bead while 
two counter-propagating fronts, coherent phonons, travel 
through the chain. We found that the energy transport 
is identical to that of a linear chain, with either mass or 
coupling disorder, at least up to the time scales reached 
in our simulations. The spreading of the energy is char¬ 
acterized by an asymptotic time dependence of the mean 
second moment of the energy m 2 of the form {m 2 ) ^ 
where = 0.5 and /3 = 1.5 for an initial displacement 
and an initial momentum excitation respectively. Ad¬ 
ditionally, the mean participation number (P) in this 
regime, converges to a constant value, which depends 
on the strength of the disorder. 

For larger values of the initial conditions, in the weakly 
nonlinear regime, we found that for initial displace¬ 
ment excitations the energy spreading does not exhibit 
a clear asymptotic time dependence. However it does 
cross to a super-diffusive behavior, since for large enough 
timescales, the slope of m 2 becomes larger than 1. In fact 
this behavior, which was also observed in the recent study 
of [26] for an FPU lattice, is found to be closely connected 
with the nonlinear dynamics of the localized state formed 
around the central spherical bead. We found that after 
a sufficient time interval (between 10 ^ and 10 ^ normal¬ 
ized time units; see Eq.Q), the central localized region 
consisting of about 10 beads, starts to delocalize and the 
energy stored in these spherical beads starts to radiate 
into the system. In this weakly nonlinear regime, the 
dynamics is governed by the power Hertzian nonlinear¬ 
ity. On the other hand, for initial momentum excitations, 
the slope of m 2 remains around 1.5, as in the near hn- 
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FIG. 11. (Color online) Profiles of the energy distribution for initial displacement (top panels) and momentum (bottom panels) 
excitation. The values of the initial displacement excitations are u — 0.01,1,10 for the first three top panels (frome left to 
right) and u — 1 when F — {) for the fourth top panel. The values of the initial momentum excitations are ii — 0.01,1,10 
for the first three top panels (from left to right) and ii — 1 when F — ^ for the fourth top panel. These values cover all the 
dynamical regimes: from the near linear (left column) to the highly nonlinear (right column) regime. The different profiles are 
taken at times t ~ 1500 (blue), t ~ 5000 (green) and t ~ 10000 (red). The dashed lines denote slopes of 5/2 (top) and 3/2 
(bottom). For the top third panel, the additional solid line denotes the slope of 3/2. 


ear regime, but (P) exhibits the same behavior as for 
displacement excitations. 

For even larger amplitudes of the initial excitation, the 
energy transfer becomes substantially different. The 
energy profile of the chain reveals an almost ballistic 
behavior with an almost equal distribution of the en¬ 
ergy around the excited portion of the chain. In this 
regime, which we characterized highly nonlinear, the sys¬ 
tem exhibits a large number of opening of gaps between 
beads, which is a nonsmooth nonlinear process. We 
found that these gaps propagate in the chain, and that 
the transport of energy is mediated by a shock-like struc¬ 
ture, which bares similarities with the selfsimilar solution 
found in [58] . 

An important result of our work is the following: al¬ 
though it is known that the energy transport for the dis¬ 
ordered linear chain strongly depends on the type of ini¬ 
tial conditions (i.e. displacement or momentum excita¬ 
tions), we found that in the highly nonlinear regime it is 
independent of the initial condition. This is a rather gen¬ 


eral feature of disordered granular chains as it was very 
recently observed in different disordered dimer granular 
chain [52]. In particular, in the linear case the asymp¬ 
totic time dependence of (m 2 ) shows a slope of P = 0.5 
(displacement) and (3 = 1.5 (momentum), while for the 
extreme nonlinear limit of P = 0 both initial conditions 
lead to a slope of (3 ~ 1.8. Additionally, the energy pro¬ 
files in this regime show that there is no distinct localized 
state. 
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